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1 Introduction 

Elasticity has the reputation of being a rather boring and un-physical sub- 
ject. In fact, dealing with second (stress-strain) or fourth (elastic constants) 
rank tensor guarantees that the notations are in general rather heavy, and 
that the underlying physics is not easily captured 1 . The elastic stress-strain 
behavior is, however, a very basic property of all solid materials, and one that 
is rather easy to obtain experimentally. Hence simulation methods for solid 
systems (hard or soft) should in general consider the obtention of these elas- 
tic properties. The determination of the response of a system to an imposed 
external stress is an important topic, which was first addressed in the sim- 
ulation community through the pioneering work of Parrinello and Rahman 
|21 El- While the main issue addressed by the Parrinello- Ray- Rahman method 
was that of allotropic transformations in crystalline solids, the current inter- 
est in nanostructured materials opens a number of new questions. What are 
the appropriate measures of stress and strains at small scale? Can one scale 
down the constitutive laws of macroscopic elasticity, and if yes to what scale 
? Are the elastic constants of nanometric solids identical to those of the bulk? 
What about the elastic/plastic transition at small scales ? Many of these 
questions are still unanswered, and the object of active studies. In the fol- 
lowing, I will describe recent studies that investigate some of these questions. 
After briefly recalling the appropriate definitions, I will particularly concen- 
trate on the case of amorphous systems, in which some surprises arise due to 
the non-affme character of deformations at small scales. In the last section, I 
will describe how the Parrinello-Rahman scheme can be extended to systems 
in which the particle representation has been replaced by a field representa- 



1 To appreciate the very nice physical content of elementary elasticity theory, the 
reader is referred to the "Feynman lectures on Physics" or in a more formal style 
to the first chapters of the Landau and Lifschitz textbookpQ 
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tion, as is common in polymer or amphophilic systems that self organize at a 
mesoscopic scale. 



2 Some definitions 

A convenient way to study elasticity with a microscopic viewpoint is to think 
of systems contained in periodic cells with variable shape. This approach, 
which was initiated by Parrinello and Rahman, is useful from a practical, but 
also conceptual viewpoint. In a system with periodic boundary conditions, the 
simulation cell can be defined by three (in general nonorthogonal) independent 
vectors hi,h2,h3 forming the sides of the parallelepiped cell. The Cartesian 
coordinates of these vectors can be used to construct a 3 x 3 matrix h defined 
by h = (hi, ti2, 113). The Cartesian coordinates of any point R in the cell can 
be expressed as 

R = ftX (1) 

where X is a rescaled vector whose components lie in [0, 1]. Integrals on R 
can be converted into integrals over X by using a scaling factor det ft,, which 
represents the volume of the cell, V. In the case of a particle or monomer 
number density, for example, one can write 

p(R) = p(X)(deth)~ 1 (2) 

The metric tensor G is constructed from h as 

G = h T h (3) 

where ft T is the transpose of h. G is used in transforming dot products from 
the original Cartesian to rescaled coordinates, according to 

R • R' = X • G • X' = X a G a fi Xa (4) 

where here and in the following summation over repeated indexes is implicit. 

Elasticity theory describes the deformation of any configuration from a 
reference configuration in terms of a strain tensor. This tensor is constructed 
by relating the vector connecting two points in the deformed configuration to 
the corresponding displacement of the same points in the reference configu- 
ration. If the reference configuration of the simulation box is denoted by h , 
the displacement is u = R — Ro = (hhn 1 — l)Ro, and the strain is given by 
citcParrincllo8 1 ,Ray84 
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{!?o)- l h T h(h )- x -l =- (^or^o)- 1 -! (5) 



where 1_ denotes the unit tensor. It is important to note that this expression, 
usually known as the Lagrangian strain tensor is not limited to small defor- 
mations pQ. Usually, the reference configuration h will be defined as a state 
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of the system under zero applied external stress. If one starts with a cubic 
cell, h Q is the identity matrix and the relation between e and G simplifies. The 
thermodynamic variable conjugate to this strain tensor, in the sense that the 
elementary work done on the system can be written in the form 

SW = V Tr(t 6e), (6) 

is the thermodynamic tension tensor t 0], also known as Piola-Kirchhoff 
second stress tensor. Vb = dcth denotes the volume of the system in the 
reference configuration. This thermodynamic tension tensor can be related to 
the more usual Cauchy stress tensor a through 

3. =yh(h o r 1 t_(h T )- 1 h T (7) 

The tension is the derivative of the free energy with respect to the strain, which 
is calculated from the reference configuration. The Cauchy stress, on the other 
hand, is the derivative of the free energy with respect to an incremental strain 
taken with respect to the actual configuration. This Cauchy stress tensor is the 
one that enters momentum conservation and whose expression is given by the 
usual Irving-Kir kwood formula for pairwise additive potentials (see below). 
The difference between these two quantities can be understood, qualitatively, 
from the fact that the strain is not an additive quantity, as can be seen from the 
existence of the nonlinear term in equation [5] While the Cauchy stress has a 
mechanical meaning in terms of forces within the sample, the thermodynamic 
tension is a purely thermodynamic quantity, and does not in general have a 
simple mechanical interpretation. 

Fortunately, in the limit of small deformations which I will concentrate 
on, the differences between these various expressions of stress tensors can be 
forgotten. This is not the case, however, for large deformations, where these 
differences result in a whole variety of stress/strain relations and associated 
elastic constants. This is especially important when dealing with solids under 
high pressure, where for example one has to be careful as to which of these 
elastic constants is used to compute e.g. sound velocities. I will refer the 
interested reader to the reference publication of Klein and Baron jS] for an in 
depth discussion of these subtleties. 



3 Finite temperature elastic constants: Born and 
fluctuation terms 

The elastic constants for a material made of particles interacting through a 
pair potential 4>{r) (I'll keep this simplifying assumption in the following) 
can be determined from simulations using an approach presented by Hoover 
and coworkers . These authors start from the explicit expression of the free 
energy in terms of a configuration integral 
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exp(-pF) = V N J dX 1 dX 2 ..dX N exp(-/3#({/iXi})) (8) 

where i?({R.J) = J2ij — R3) i s the total interaction energy. 
The derivative with respect to strain is taken using 

BF f BF \ 

dF = Tr(—dG)=2T±\h —h T de) (9) 

which gives for the thermodynamic tension matrix 

V t a/ 3 = Nk B Th a _ ay G~gh T 'o y sp+ho,arr(^2 x i3,f^ij,s - L ^ ))ka,80 = Nk BTh a7 G^gh T Sl3 +(f a p) 

*3 

(10) 

where = Xi —Xj, yV 1S the summation over all distinct pairs of particles, 
and the pair potential <f> is assumed to depend only on the particle separation 
Rij. The brackets () denote a thermal average. This obviously reduces to the 
usual Kirkwood formula for small deformations (ft. = h). Note that the first 
term arises here from the volume factor V N in the configuration integral, but 
could also be obtained by introducing the momenta and the kinetic energy 
contribution in the partition function. The last term defines the potential en- 
ergy contribution to the microscopic stress tensor, denoted by T_. Carrying out 
one more derivation with respect to strain, one obtains the elastic constants 
in the limit of zero strain (more general expressions for arbitrary strain can 
be found in refs 011]): 



Ca^s — Q a/3 = 2Nk B T(5 ai 5ps+5 a s5i} 1 )-— %- (T a pf 7 g) — (T aj g}(T 7 i) 
Ce 7 a kb-L l 

(11) 

Where the last (Born) term is written in terms of potential energy functions 

si B om 

ij \ l J l 3 / 

with Rij = Ri — Rj . The term in square brackets in ^2 is called the fluc- 
tuation term, and is generally expected to be a correction to the main Born 
term at finite temperature. We will see below that this term, even in the low 
temperature limit, remains an essential contribution to the elastic properties 
of disordered systems. Formulae that generalize the equations above to local 
stresses and elastic constants can be found in references [3] E| • 

Note that a different approach must be used for hard core potentials. In 
that case, one possibility is to study strain fluctuations either in a cell of 
variable cell m, or to directly study the stress-strain relation in a deformed 
cell H2|. 



s~iBorn 
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4 Amorphous systems at zero temperature: nonaffine 
deformation 

The simulation of amorphous systems at low temperatures is interesting from 
the point of view of the physics of glasses, but also because these models can 
serve as very elementary examples of "athermal" systems such as granular 
piles or foams. The elastic/plastic response of such complex systems is not 
well understood, and is currently the subject of many experimental studies 

[ni m nuns]. 

A naive approach to the calculation of elastic properties in such systems 
would consist in taking the second derivative of the potential energy H = 
J2i<j 'fi(Rij) with respect to strain. Such an approach is easily shown to yield 
elastic constants that correspond to the Born expression, without the thermal 
average brackets. Although such an approach seems natural - the "fluctuation" 
term could be ignored at zero temperature -, it proves in fact completely 
incorrect for disordered systems, or even for crystals with a complex unit cell. 

The essential point is that the derivatives have to be taken not at constant 
X, but rather keeping the force on each atom equal to zero in the deformed 
configuration. In other words, one has to allow for relaxation of the deformed 
configuration before computing the energy and stresses. It was shown by Lut- 
sko [7] (see also the recent work by Lematre and Maloney, ref|5]) that this 
relaxation gives a contribution to elasticity which is identical to the zero tem- 
perature limit of the fluctuation term. The corresponding proof can be briefly 
summarized as follows. 

The elastic constant is written as 

C a f3yS = q \F i= (13) 

pFi = indicates the constraint that forces on the particles must remain 
zero during the deformation. The variables in the problem are the reduced 
coordinates, X^, and the strain e . The force F.; is a function of these variables, 
so that the constrained derivative above can be written as (for simplicity, 
we drop in this formula and the following the Greek indexes for Cartesian 
coordinates): 

dt dt ( 8Fj \ _1 dFj 

Q= ^-m{dxJ ~dT (14) 

where terms such as f§3c J nave to be understood in a matrix sense. In fact, 

Dij = (fit 2 ") i s nothing but the dynamical matrix of second derivatives of the 
potential energy with respect to atomic positions. Finally, using the definition 
of the force Fj and of the tension t as derivatives of the potential energy, 
equation 1141 can be rewritten in the more symmetric form 
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The direct evaluation of the second term in this equation is not straight- 
forward, hence the actual procedure to obtain the zero temperature elastic 
constants generally consists in carrying out explicitly an affine deformation 
of all coordinates, then letting the atomic positions relax (using e.g. conju- 
gate gradient minimization) |22J to the nearest energy minimum. Eauationll5l 
however can be used to show that the resulting elastic constants are identi- 
cal to those obtain at a finite, low temperature using equation II II The proof 
goes simply by expanding both the stress and energy in terms of the atomic 
displacements in the unstrained reference configuration 

1 dt 
H = H + -D lJ SX l SX J ■ t = t + —SX l (16) 

Performing the resulting gaussian integrals, it is seen that the "fluctuation" 
term in 1111 and the "relaxation" term in ^1 are identical in the limit of zero 
temperature. 



5 Numerical results 

We now address the importance, qualitative and quantitative, of this "relaxation- 
fluctuation" contribution. Quantitatively, the importance of this contribution 
can be estimated from figure^ hi this figure, the Lame coefficients of a three 
dimensional, amorphous, Lcnnard- Jones system (see ref. |17p at zero tem- 
perature are computed using the Born approximation and the exact formula, 
eauation llll It is seen that the relaxation term can account for as much as 50% 
of the absolute value of elastic constants. Although this fraction may obvi- 
ously be system dependent, the situation is very different compared to simple 
crystals (with one atom per unit cell, e.g. FCC in the Lennard- Jones system) 
in which the elastic constants are exactly given by the Born term (see e.g. |18| 
for a comparison between amorphous systems and simple crystal structures) . 
The relaxation contribution tends to lower the shear modulus /i, and to in- 
crease the coefficient A. Remarkably, the bulk modulus K = A + 2/i/d w 57 
(d = 3 is the dimensionality of space) would be correctly predicted by the 
Born calculation. 

As discussed above, the Born formulae would be exact at zero temperature 
if the global deformation was equivalent to an affine deformation of atomic 
coordinates at all scales, i.e. a mere rescaling of h at fixed values of {Xi}. The 
failure of the Born calculation can therefore be traced back to the existence 
of a non-affine deformation field, which stores part of the elastic deformation 
energy. This field is defined by substraction from the actual displacement 
of the atoms (after relaxation) the displacement that would be obtained in 
the affine hypothesis. The existence of this non-affine deformation field was 
pointed out in several recent publications E3 1201 HH hi Ell H3- it was 
in particular shown that this non-affine contribution is correlated over large 
distances, and is organized in vortex like structures (i.e. is mostly rotational 
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Fig. 1. Lame coefficients A (spheres) and fi (squares) vs. system size L, for a sim- 
ple polydisperse Lennard- Jones "glass". Full symbols correspond to a direct mea- 
surement using Hooke's law with relaxation, open symbols correspond to the Born 
approximation). The effect of system size is weak. For large boxes we get m 15 
and A « 47. 



in nature). These properties are illustrated in figures [21 and for a simple 
Lennard- Jones two dimensional system. 

A slightly different, more local and general, definition of the nonaffine dis- 
placement field (or "displacement fluctuation") was proposed in ref. [5]. In 
this reference, the nonaffine field is defined by substracting from the actual 
displacement a local displacement field built that is obtained using a coarse- 
graining procedure. This allows in principle to deal with situations in which 
the displacement field has a complex structure. In the case of simple shear 
considered here, our definition should be sufficient. Reference [5] also demon- 
strates that, even when the displacements are not locally affine, there exists 
a local linear relation between the stress and strain fields, at sufficiently large 
scales (resolution) . These fields are evaluated at the same position with a cho- 
sen resolution. The derivation does assume that the displacement fluctuations 
are uncorrelated over sufficiently large scales, i.e. it will be valid at scales 
larger than the one discussed here. 

A quick study of a simple one dimensional model is useful to understand 
the importance of the nonaffine displacement field. Let us consider a chain 
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Fig. 2. Snapshot of the nonaffine displacement field in a 2d Lennard- Jones amor- 
phous system undergoing uniaxial extension. Note the large scale, vortex like struc- 
tures. The sample contains about 20 000 particles. 

of N atoms, connected by springs ki, submitted to a force F. The extension 
of spring i (linking site i to i + 1) is Si = F/ki. One can therefore write the 
displacement of atom p 

u^Fxj^ki 1 (17) 

i=l 

The affine displacement is just u a J s = (p/N) xFx Yn=i K 1 = P F < k ^ 
where the <> refer to an average over the distribution of elastic constants 
and the large N limit has been taken to compute the affine displacement. As 
a result we have for the nonaffine displacement of atom p 



Microscopic elasticity of complex systems 



9 



0,02 



O 0,00 

~¥ 
V 

^-0,02- 

o 

3 -0,04 



V "0.06 



-0,08 



o — e 2D : L=104 
3D : L=64 



W 3D A 2D 



Anti-Correlation ! 



_L 



_L 



10 



20 



30 40 

r/c 



50 



60 




Fig. 3. Correlation function C(r) =< u 



(r) > of the nonamne displace- 



ment field u" A (r) in an amorphous sample undergoing a simple uniaxial extension. 



Both the 2d and 3d case show correlations that extend over scales of typically 20-30 
particle sizes. A negative tail can be associated with "vortex like" structures that 
reflect the essentially rotational character of the non affine field. 



P 

u% A = u p - uff —Fx ^(fcr 1 - < k- 1 >) (18) 

i=l 

which shows that in this simple Id situation the mean squared value <( 
Up A ) 2 >oc p < (Sk^ 1 ) 2 > is increasing linearly with p, and proportional 
to the variance of 1/ki (see also the discussion by DiDonna and Lubensky, 
cond-mat/0506456). 

The existence and nature of the length scale over which the non affine field 
is correlated is still a matter of debate 2 . Clearly, the correlation length £ is 
a lower limit for the applicability of continuous elasticity theory. This limit 

2 In a recent preprint ( cond-mat/0506456 , "Nonamne correlations in Random Elas- 
tic Media"), DiDonna and Lubensky argued that the nonamne field has logarith- 
mic (in 2d) or 1/r (in 3d) correlations, and hence no characteristic length scale. 
That such singular behaviour is possible is already illustrated in the simple Id ex- 
ample above. Such a behaviour, however, is not evident in our numerical results. 
Their calculation, based on the fact that elastic propagator have 1/fc 2 behaviour 
in Fourier space, is perturbative in the disorder strength, and it could be that we 
are investigating a "strong disorder" limit. In any case, further investigations are 
needed to assess the actual existence of such long range correlations 
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manifests itself in several different ways. If one considers the vibrations of a 
system of size L, these vibrations will be properly described by the classical 
elasticity theory only if the corresponding wavelength is larger than £. If one 
considers the response to a point force, this response will be correctly described 
by the continuum theory only beyond the length scale £. More precisely, it was 
found that the average response is described by continuum theory essentially 
down to atomic size, but that the fluctuations (from sample to sample) around 
this average are dominant below £ [22] . 

Finally, it is very likely that the existence of this length scale is related to 
a prominent feature of many disordered systems, the so called 'boson peak'. 
This feature actually corresponds to an excess (as compared to the standard 
Debye prediction, g(oj) cx in d dimensions) in the vibrational density of 
states g(u>) of many amorphous systems. This excess shows up as a peak in a 
plot of <7(o;)/w d_1 vs u, that usually lies in the THz range. In terms of length 
scales, we found that this peak typically corresponds to wavelengths of the 
order of magnitude of £ (see figure 0J. A simple description [23] is therefore to 
assume that waves around this wavelength are scattered by inhomogeneities, 
and see their frequencies shifted to higher values. Pressure studies show that 
the boson peak is shifted to higher frequencies under pressure, consistent with 
a shift to smaller values for £ obtained in simulations. Another very interesting 
evidence for the existence of mesoscale inhomogeneities was recently provided 
by Masciovecchio and coworkers [2Jj , by studying Brillouin spectra in the ul- 
traviolet range. The width of the Brillouin peak shows a marked change for 
wavelength between 50 and 80 nm, indicative of scattering by elastic inhomo- 
geneities. 

A very interesting question is whether this characteristic correlation length 
for elastic inhomogeneities, which can reach rather large values compared to 
atomic sizes, is somehow associated with a 'critical' phenomenon. An idea that 
was recently suggested by Nagel and co-workers is that this correlation 
length should diverge at the so-called 'jamming' transition in purely athermal 
systems. The jamming density is defined, in a system with purely repulsive 
interactions at zero temperature, by the density at which the system will start 
to exhibit mechanical rigidity. Below the jamming density <p c , an infinitesimal 
temperature results in diffusion, while systems above </> c remain in a frozen 
state on macroscopic time scales. Based on general arguments concerning iso- 
staticity of the packing at <f> c , Nagel and coworkers suggested the existence 
of a correlation length associated with soft modes, that diverges at the tran- 
sition. Although the arguments are in principle valid only for contact type 
interactions, it would be quite interesting to follow the evolution of £ for a 
system with attractive interactions, but under tension, expecting perhaps a 
divergence close to the rupture threshold. 

Finally, let us mention that a different way of studying local elastic prop- 
erties was proposed by de Pablo and coworkers through the study of local 
elastic moduli, which can be defined by using the definition ^2 to a small, 
finite box ^Hj- Depending on the scale at which they are measured, these 
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Fig. 4. Vibrational density ol states in a 3d Lennard- Jones amorphous system. The 
"boson peak" is apparent as the deviations from the Debye prediction in the ratio 
g(uj)/u> 2 . This peak is observed at frequencies ol order 2nc/£, where c is the speed 
of sound. 



moduli can take negative values. Such regions would be unstable, if they were 
not immersed in a matrix of " normal" regions. The size over which local elas- 
tic constants are found to be negative is small (typically 3 particle sizes) , and 
could probably be considered as a first coarse-graining scale for using classical 
methods for disordered systems 



6 Polymeric systems: stresses and self consistent field 
theory 

Self consistent field theory is a powerful approach to the determination of 
phase equilibria in polymer systems with complex architectures. The theory 
directly deals with density fields rather than particles, and minimizes a mean 
field like free energy. The method is particularly suitable for polymers, in 
which interactions on length scales comparable to the chain size are effec- 
tively very soft. The method is well known and has been described in many 
publications (see e.g. |211|2SI) and I will only describe briefly the main steps 
and the way a constant stress method can be introduced in the simulation. 
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This option allows one to obtain relaxed configurations at zero imposed stress 
easily, or to study the effect of anisotropic tension on phase behavior. 

As a representative example, I consider a model for an incompressible 
AB diblock copolymer melt |2§]- The melt consists of n identical diblock 
copolymer chains composed of monomer species A and B and is contained 
in a volume V . Each of the chains has a total of N statistical segments; a 
fraction / of these segments are type A and constitute the A block of each 
macromolecule. For simplicity, the volume occupied by each segment, vo, and 
the statistical segment length, b, are assumed to be the same for the A and B 
type segments. The Hamiltonian for this system can be written 

H = ib^J o lds (^r) 2 + v 0X ABk B T I dvp A {v)p B {v) (19) 

where Rj(s) with s S [0, 1] is a space curve describing the conformation of the 
ith copolymer and R 2 g0 = b 2 N/6 is the radius of gyration of an ideal chain of 
N statistical segments. Interactions between dissimilar segments A and B are 
described by the Flory parameter xab- The densities pa,b{ v ) are microscopic 
segment density fields defined by 

" rf 

p A (r)=AV/ ds<5(r-R,( S )) (20) 
i=i Jo 

and 

n -i 

p B {v) = AV / ds 6(r - Ri(fl)) (21) 
<=i J f 

A local incompressibility constraint /5a (r) + Pb{?) = Po is imposed in this 
standard copolymer melt model for all points r in the simulation domain. The 
total segment density po can evidently be expressed as po = nN/V = 1/vq. 
Using the rescaled coordinates X(s) (taken in [0, l] 3 ), the generalized partition 
function that has to be sampled for a fixed value of the thermodynamic tension 
t reads 

Z = J d(h){deth) nN S{deth-V Q ) exp{-f3V t : e) 

n - 

x II / 2 ? X l ( S )exp(-/3i/)[] ( 5( / 5 A (x)+p B (x) - nN) (22) 

(t : e is the contraction t a pea(3). The final factor in the above expression 
imposes the constraint of local incompressibility. Moreover, incompressibility 
implies globally that the cell volume remains fixed at its initial value, i.e. 
deth = V = Vo = deth . This is enforced by the delta function in the first 
line above. Hence all shape transformations should be volume preserving. The 
practical implementation of this constraint will be discussed below. 
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Hubbard-Stratonovich transformations are used to convert the particle- 
based partition function [221 into a field theory [2B] • These can be carried out 
straightforwardly on the polymer partition function for a given cell shape h, 
Z(h), with the result 

Z(h) = l[ / DX,(«)exp(-/3ff)IJ<J(p A (x) + p B (x) - nN) 

i=l J x 

= J Vw exp(n In Q[w,K\- E[w]) (23) 

where Q[w,h] is the partition function of a single copolymer chain experienc- 
ing a chemical potential field u>(x, s), J Vw denotes a functional integral over 
the field w, and E[w] is a local quadratic functional of w that reflects the A-B 
monomer interactions and the local incompressibility constraint: |28j 



71 

E[w] = - / '/x 



(w B - w A ) 2 ~ {w A + W B ) 



2 X N 



(24) 



Here we have noted that for an AB diblock copolymer melt, the potential 
w(x, s) amounts to a two-component potential, i.e. w(x.,s) = wa(x) for s G 
[0, /] and u>(x, s) — iub(x) for s e [/, 1]. 

The object Q[w, h] is a normalized partition function for a single copolymer 
experiencing a potential field w(x, s) This partition function can be obtained 
from a single-chain propagator g(x, s) that is the solution of a modified diffu- 
sion equation 

subject to g(x, 0) = 1. The single chain partition function is given by Q[w, K\ = 
Jdxq(x, 1). 

Finally, the partition function for an incompressible diblock copolymer 
melt confined to a cell of variable shape can be expressed as a field theory in 
the variables h and w: 

Z = J d(h)(deth) nN J Vw S(deth- V ) exp(-F[w,h\) (26) 

where F [w,h\ is an effective Hamiltonian given by 

F[w,h]=/3V t :e + E[w]-n]nQ[w,h\ (27) 

In the mean- field approximation (SCFT) , for a given shape h of the simulation 
box, we approximate the functional integral over w in eq 1261 by the saddle 
point method. For this purpose, the functional Q[w,h\ can be evaluated for 
any w and h by solving the modified diffusion equation (using e.g. a pseudo- 
spectral approach). The saddle point (mean-field) value of w, w* , is obtained 
by applying a relaxation algorithm [2H1 00] to solve 
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SF[w,h] 
5w(x, s) 







(28) 



In the mean-field approximation, F[w*,h\ corresponds to the free energy of 
the copolymer melt (in units of kgT). 

In a simulation at constant tension, the relaxation equation for the fields 
must be supplemented by a corresponding evolution equation for the cell. This 
equation is chosen to be a simple relaxation 



dh 
~dt 



-A hDh 



-idF[w,h] 
dh 



(29) 



where the tensor D_ is a projection operator whose action on an arbitrary 
tensor Mis a traceless tensor, i.e. D_M_ = M_ - (l/3)Tr(M )1. Equation 1391 
corresponds to a cell shape relaxation that (for Ao > 0) is down the gradient 
dF/dh, approaching a local minimum of the mean-field free energy F[w* ,h\. 
The "mobility" tensor hDh 1 is chosen so that the cell shape dynamics de- 
scribed by ea l2*9l conserves the cell volume. 

Application of ea 1291 requires an expression for the thermodynamic force 
dF I dh. Explicit differentiation, noting the constraint of constant det h, leads 
to 

dF[w,h] _ / d 



dh 



PV 



-Tr^+MZ 



(30) 



where £_ is a symmetric tensor defined by 



E a p[w,h] 



2k B TndlnQ{w,h] 
V 

ksTn 



2VR% 



dG a p 

1 dfj dX a (s) dXp(s) 







ds 



ds 



(31) 



The angular brackets in the second expression denote an average over all 
conformations X(s) of a single copolymer chain that is subject to a prescribed 
chemical potential field w and fixed cell shape h. 

The first term on the right hand side of eo l30l can be conveniently rewritten 

as 



)) = ^-(Tr V l h T hh^H ) = ihh^H h T \ p . 
Hence, eq|5!5]can be compactly expressed as 



(32) 



dt 



{h^t h 1 



o )+£ 



(33) 



where A > is a new relaxation parameter defined by A = (3Vq\q. 

Equation 1331 will evolve the cell shape to a configuration of minimum free 
energy (in the mean-field approximation). This configuration can either be 
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metastable (local minimum) or stable (global minimum). Addition of a noise 
source to the equation provides a means for overcoming free energy barri- 
ers between metastable and stable states, i.e. a simple simulated annealing 
procedure. 

An equilibrium solution of the cell shape eauation l33l is evidently obtained 
when 

ihohh^o 1 ) + S =0 (34) 

Combining eqs 1311 and 1341 it is seen that this equilibrium condition corre- 
sponds to a balance between the externally applied Cauchy stress, a , and the 
internal elastic stress, a , sustained by the polymer chains 



(35) 



where 



This expression for the internal polymer stress is well-known in the polymer 
literature [3*T] . 

Eauation l33l drives a change in the shape of the simulation cell (at constant 
cell volume) to approach the equilibrium condition 1341 at which the internal 
elastic stress of the copolymers balances the imposed external stress. 

The last step is to find an expression for the internal stress tensor g_ mt 
(eq 13611 or S_ (eq 13111 in terms of the single chain propagator, which is the 
central object computed in a field-theoretic simulation j^H]. The appropriate 
expression turns out to be 

■» - ^f^/' ^Wii) .. (37) 



(n/V)k B T Q ia J J dX 7 dX s &(i 

we refer to the appendix and to |32j for the derivation of this expression. 

From the above derivations, it is clear that the quantity that is exter- 
nally imposed in the method is not the Cauchy stress, but rather the thermo- 
dynamic tension. The Cauchy stress, which is the experimentally accessible 
quantity, is a result of the simulation, as is the cell shape. This feature is 
general in any application of the Parrinello- Rahman method in which a par- 
tition function of the form refzpol is sampled using Monte-Carlo, Langevin 
or molecular dynamics. The Cauchy stress has therefore to be obtained inde- 
pendently, using equation ^] in the present case, or in the case of molecular 
system through the Irving-Kirkwood formula. 

To illustrate the method, figureElshows the evolution of the simulation cell 
under zero tension in a simple case. In the system under study, the parame- 
ters have been chosen so that the equilibrium phase is ordered on a triangular 
lattice, under zero external stress. Starting with a square simulation cell, evo- 
lution to the correct rhombohedral shape is obtained after a few relaxation 
steps. Other examples of application may be found in |32| . 
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Fig. 5. Transformation of a square cell under zero external stress, when the melt is 
quenched into the stability region of the cylindrical phase (\N — 15.9, / = 0.64). 
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A Expression for the stress tensor in SCFT 

To obtain equation 1371 we start with the definition of J7 

m/r 9 1 dQ[w,G] 

PVE «e = Q~dG^T~ (38) 

The derivative of the single chain partition function can be calculated by 
discretizing the paths with a small contour step A. 
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2 dQ[w,G] 1 pi , r , Y f , Y , 

J 0X(a)*(X - X(a))*(X' - X(,s + Z\)) (^3^) 
«P Jo da^Gap^. - £ ds W (X(s), a))39) 

Except between the points X, s and X', s+Z\ one can replace the path integrals 
with propagators q, so that 



2 dQ[w,G] 

( exp (-^Getf (X a - - X' p ) 



Q ~dG^T = ^oQ & dsJdXl dX ' q{X > s)q{X ' ' 1 ~ s ~ A) 



One can then set X' = X + u, and expand for small u and A according to 
9 (X + u, 1 - s - A) = 9 (X, 1 - -) - A d J- g + + (41) 

The derivative w.r.t. s can be eliminated by applying the modified diffusion eq 
1251 One also requires second and fourth moments of the Gaussian distribution 
of displacements u, 

u^up = G-'(2R 2 g0 A) 
u a u p u lU& = {2R 2 g0 Af{G-lG-j + G-]G- l + GjG^) 
By means of these results, we have 



^ I ^ + G-/dX w (X)p(X) 



2 dQ[w,G\ _ 
Q 9G a f3 

+ %G- a lG-l + G-} t G-l + G-jG-i) J dxfi dsq(X, s)^ 



(40) 



dxs 

(42) 



where p(X) = Q 1 f Q ds g(X, s)q(K, 1 — s) is the single-chain total monomer 
density operator. There is a partial cancellation in the last two terms so that 

2 dQ[w,G] i/dxppc) j , Y , 

= G Q/3 ^^ + iGQ/3 J dXw (XV(X) 

+ ^^G^/dX^^X,,)^^ 

(43) 



The internal polymer stress is obtained after matrix multiplication by h 
on the left and ft. T on the right. This implies that the first two terms be- 
come a simple isotropic stress contribution, and are therefore not relevant to 
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an incompressible system. The final formula for the internal stress tensor is 
therefore, apart from this diagonal contribution, 



The tensor S appearing in equation 1331 is given by an expression similar to 
1441 with G replacing h. The factor ksTn/V accounts for the total number 
of chains, and produces a stress with the correct dimensions. In practice, 
the stress will be made dimensionless by dividing by this factor, so that the 
dimensionless stress is given by 



Equation 1371 is obtained after integrating by parts. 

A local (rather than volume averaged) version of this connection between 
the stress tensor and the polymer propagator was derived previously in |33j . 
Numerically, <j a p is evaluated from ea 1451 using a pseudo-spectral scheme. The 
derivatives with respect to spatial coordinates are obtained by multiplying the 
propagator by the appropriate components of the wavevector in Fourier space, 
and transforming back into real space. 
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